ポリゴン内にある点の集計

sp パッケージの over 関数を用いてポリゴンの内点を判定し、内定んデータで各種統計値を算出する。
例として、統計局の 地図で見る統計（統計ＧＩＳ） からダウンロードした広島市中区の世界測地系の地図データに対して、
乱数で点データを生成し、各区域ポリゴンの内点の個数を集計し可視化する。

library(maptools)

# 世界測地系
pj <- CRS("+proj=longlat +datum=WGS84") 

# ダウンロードしたshapeファイルのリスト
d <- readShapeSpatial("map/h22ka34101.shp",proj4string=pj)  # 町丁・字境界データ
d <- subset(d, d@data$HCODE == 8101)  # 水上の境界を除外

# 地図を描く
plot(d)

# 乱数で点データを発生
n <- 100
pnt <- data.frame(
id=1:n,
x=runif(n=n,min=d@bbox[1,1],max=d@bbox[1,2]),
y=runif(n=n,min=d@bbox[2,1],max=d@bbox[2,2])
)

# SpatialPointsDataFrameに変換
coordinates(pnt) <- ~x+y
pnt@proj4string <- pj
points(pnt,pch=19,cex=0.2)

# 各ポリゴンの内点の個数
d@data$n.inner <- over(d,pnt,fn=length)$id
d@data$n.inner[is.na(d@data$n.inner)] <- 0
table(d@data$n.inner)

## 
##  0  1  2  6 
## 86 22  3  1

# 内点の個数で色分け地図を作成
spplot(d,zcol="n.inner")